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ABSTRACT 


This thesis deals with the study and some applications 
of the Method of local Variations which is a nei^bourhood 
method utilising a direct search technique for obtaining 
locally optimal solutions to optimal control problems. A 
brief review of direct and indirect methods for solving 
optimal control problems has been given, together with an 
introduction to the Method of Local Variations technique 
and the motivation for its study. This is followed by a 
detailed description of the local variations algorithm, 
its comparison with dynamic programming techniques, and 
some of its advantages and disadvantages, Next, the 
Method of Local Variations and some modifications on it 
have been applied to solve several optimal control problems. 
Consequently, some important features of the local variations 
algorithm have emerged, providing a basis for its comparison 
with the other better known techniques of optimisation and 
in drawing conclusions about its relative merits and demerits. 
Finally, the conclusions arrived at, together with some 
suggestions regarding modifications to this method for 
improvement of the basic algorithm, have been given. 
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CHAPTER 1 
IHTRODHCTION 

A problem which seeks to extremis e a function 
(or functional) of one or more variables and in which 
the solution is required to satisfy certain specified 
constraints, may be referred to as an optimisation 
problem. A d 3 mamical system may have such, a problem 
defined for it; in the field of control theory, this 
problem is known as an optimal control problem. 

The formulation of an optimal control problem 
requires: 

a) a mathematical model of the process to be 
controlled, 

b) a statement of the physical constraints and 

c) the specification of a performance criterion. 

Attention will be confined, in this study, to dynamical 
systems which can be described by state equations of 
the form: 

x(t) *= * u(t),t) (l) 

where x{t) is the n-dimensional state vector, 

u(t) is the m-dimensional control vector and 
a is an n-dimensional vector function which is, 
in general, nonlinear and time-varying, with 
arguments as indicated. 



Thus, the systems considered are, in general, finite- 
dimensional, deterministic, continuous-time and 
nonlinear-time-varying. 

Considering an optimal control problem defined 
in the finite time interval [t^,t^] where t^ is the 
initial time and t^ the final time, the constraints 
may consist of the specified boundary conditions on 
the state variables together with any state and control 
constraints that may be imposed in the interval [tQ,t^]. 

The performance index is a mathematical 
expression which, when extremised, indicates that the 
system is performing in the most desirable manner. Its 
choice depends on the type of problem being considered. 

¥e will consider performance indices which are expressible 
in the form: 

^f 

J * h(x(t^),t^) -f J* g(x(t), u{t), t)dt (2) 

"^0 

where h’ and g are scalar functions of the argxunents 
indicated in eqtiation (2). 

A typical optimal control problem can, thus, be 
formulated as follows: It is required to find an 
admissible control u** which causes the system described 
by equation (1) to follow an admissible’ trajectory x* 
that minimises the performance index given by equation{2), 
subject to the constraints stipulated, u* called an 
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optimal control and an optimal trajectory. The 
problems considered in this thesis all belong to the 
above category. 

There are several methods of solution of such 
problems. These methods can be broadly classified into 
two main categories, namely, direct methods and indirect 
methods. 

1.1 Direct and Indirect Methods for Solving Optimal 
Control Problems 

Indirect methods usually utilise the necessary 
conditions for optimality derived from the Calculus of 
Variations or Pontryagin’s Minimum Principle. Solution 
of optimisation problems by indirect methods can be 
obtained in the following way: 

a) a set of necessary conditions that the optimum 
must satisfy are derived; 

b) the equations constituting the necessary conditions 
are then solved to obtain candidates for the optimum 
and 

c) the candidates for the optimum are tested by using 
the necessary- and-sufficient-condition tests. 

A feature of indirect’ methods is that state and control 
constraints increase the complexity involved in obtaining 
optimal solutions. 

Application of the Cal>. — 

Pontryagin's Minimum Principle to optimal control- problems 



gives rise to’ Two Point Boundary Value Problems’ (TPBVP), 
and the difficulties associated with obtaining solution 
in this manner are due to the difficulty involved in 
solving complex TPBVP’s defined for simultaneous, 
nonlinear-time-varying differential equations. This 
aspect has been discussed, for example, by Athans [ij, 
Kirk [2], and Sage [3]. 

There are many techniques for solving TPBVP's 
encoxintered in the above formulation. The quasilineari- 
sation technique is one such method for finding a 
solution to these TPBVP’s. . Athans [l] has given a 
brief description of this technique while a more, 
elaborate treatment can be found in Sage [33* The 
quasilinearisation technique has been applied by 
Bellman and Kalaba [5, 6, 7, 8], McGill and Kenneth [9], 

Kopp and McGill [lO], and Moyer and Pinkham [ll] for 
solving TPBVP’s to obtain optimal control. This tech- 
nique has also been applied by Kumar and Sridhar [l2], 
Detchmendy and Sridhar [l3]» and Sage and Eisenberg [14] 
to solve several engineering problems in optimal control 
theory. Application of this method to discrete systems 
has been investigated by Henrici [l5]» Sylvester and 
Meyer [16], Sage and Burt [17], and Sage and Smith [I8]. 

The invariant imbedding procedure is another 
method for solving TPBVP’s. It is a technique whereby 
the missing initial (or terminal) conditions are 
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directly obtained'. An extensive treatment of this ■ 
method appears' in Sage [ 3 J while Bellman et al.[l9]» 
Detchmendy and Sridhar [20], Sage and Masters [21,22], 
and Sage and Ellis [23,24] have applied this technique- 
to solve several optimal control problems. 

« 

■ There are many other methods for solving these 
TPBVP’s. Some well known methods are the Steepest 
Descent Method, the Method of Variation of Extremals, 
the Conjugate Gradient Method, and the Second Variation 
Method. Kirk [ 2 ] has given a detailed description of 
the Steepest Descent Method and the Method of Variation 
of Extrema,ls. lasdon et al. [25] have described the 
Conjugate Gradient Method while a discussion on Second 
Variation Methods appears in Bryson and Ho [26]. 

An important feature of indirect methods is that 
a solution obtained, utilising the necessary conditions 
for optimality, may not be an optimal solution at all. 
Even if it is optimal, since the necessary conditions 
obtained from the Calculus of Variations or Pontryagin's- 
Minimxim Principle are for local optima only, the result 
will, in general, hold only locally. 

Direct methods are those which do not utilise the 
necessary conditions for optimality to obtain a solution 
to the optimal control problem. Pierre [ 4 ] defines 
direct methods of solution as those which require the 
direct use of the performance index and the constraint 
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equations of a given problem, while systematic recursive 
methods are employed to obtain an optimal solution. An 
important feature of many direct methods is that an 
initial feasible trajectory foims the starting point in 
the algorithm which finally converges to the desired 
optimal solution. At any iteration stage, the nominal 
trajectory is a feasible one. Thus, the formiilation helps 
us in picking out a good suboptimal trajectory with less 
effort than is required to find an optimal trajectory. 
Moreover, since most of the direct methods are search 
techniques, any state or control inequality constraint 
that may be specified, helps in narrowing down the domain 
of search and eases the finding of , a solution to the 
problem being considered. These methods were primarily 
devised to circumvent the problems associated with 
obtaining solutions to TPBVP's. 

An important example of a direct method is Bellman’s 
Dynamic Programming technique [27,28, 29] • This is a 
sophisticated search technique for solving multi-stage 
decision making problems. This technique is, in fact, 
the repeated, sequential, stage— by«*stage application of the 
optimality principle of Bellman. It can be applied to 
continuous time processes as well as to discrete time ones. 
Although this technique provides a global optimal control 
law, it suffers from a serious malady known as the ’curse 
of dimensionality' in that the memory requirement is very 
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large for any "but very low order problems. To overcome 
this serious drawback, several modifications have been 
made on the original algorithm giving rise to such 
techniques as Incremental Dynamic Programming of 
Bernholtz and Graham [30], Differential Dynamic Program- 
ming of Jacobson and Mayne [3l]» and State Increment 
Dynamic Programming of Larson [32]. All these techniques 
achieve a reduction in the memory requirment at the 
cost of losing the global it y of the optimal control law. 

Another example of a direct method for solving 
optimal control problems is the Gradient Projection 
Method of Rosen [ 33 » 34,55]* It is an iterative numerical 
procedure for finding an extremum of a function of 
several variables that are required to satisfy various 
constraints. A detailed description of this method in 
its application to optimal control problems appears in 
Kirk [2]. 

1.2 Motivation for the Study of the Method of Local 
Variations 

Incremental Dynamic Programming of Bernholtz and 
Graham [30], Differential Dynamic Programming of Jacobson 
and Majme [31] and State Increment Dynamic Programming of 
Larson [32] are all modified versions of Bellman’s 
Dynamic Programming technique and come \mder the category 
of neighbourhood techniques using direct methods of 
optimisation. All these techniques have one aim in 
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ccmmon - a reduction in computer memory requirement. 

•They achieve this at the cost- of losing the glohality 
of the optimal control law that they provide. 

In the incremental dynamic programming technique 
as applied to a single state variable case, a neighbour- 
hood is defined about the nominal state trajectory being 
considered. This trajectory consists of a discrete set 
of state points at certain instants of time. The dynamic 
programming algorithm is applied, in this neighbourhood 
only, to these state points, starting from the last state 
and sweeping backwards, such that the best trajectory, in 
that neighbourhood, is generated. This newly generated 
trajectory now forms the current nominal trajectory around 
which a neighbourhood is again defined and the whole 
process is repeated. This goes on till a. locally optimal 
trajectory is obtained, Por the vector case, the 
neighbourhood is defined considering perturbations in one 
state variable at a time, and the best trajectory, in. that 
neighbourhood, is generated for perturbations in that 
state variable only, keeping all other state variables 
the same. Once this trajectory has been generated for 
that particular state variable, the next state variable 
is then similarly treated. This process is repeated 
till all the state variables have been considered. This 
whole procedure is repeated till satisfactory convergence 
to an optimal control is obtained. 
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Differential dynamic programming is a successive 
approximation technique based on Bellman’s dynamic 
programming method. This technique is used for deter- 
mining optimal control for nonlinear systems. It uses 
a linearised version of the Hamilton-Jacobi-Bellman 
Equation (HJBE) at each of the time instants into which 
the time interval^ for the problem being considered, 
is divided. The optimal ■ control satisfies the HJBE. It 
is assumed that the optimal control is unknown in the 
time interval being considered, but that a nominal 
control is available in the same time period. On 
application of this nominal control a nominal state 
trajectory, in the same time interval, is produced by 
equation (l). The nominal performance index value can 
be calculated from equation (2). It is assumed that 
the optimal performance index expression is sufficiently 
well-behaved to allow a power series expansion in a 
small neighbourhood around the nominal state trajectory. 
In each iteration, the system equations are integrated 
in forward time, using the current nominal control. The 
accessory equations, which yield the coefficients of a 
linear or quadratic expansion of the performance index 
expression in the neighbourhood of the nominal state 
trajectory, are integrated in reverse time, thus 
yielding an improved control law because the performance 
index, which is a function of both state and control 
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Differential dynamic programming is a successive 
approximation technique based on Bellman's d 3 nD.amic 
programming method. This technique is used for deter- 
mining optimal control for nonlinear systems. It uses 
a linearised version of the Hamilton-Jacobi-Bellman 
Equation (HJBE) at each of the time instants into which, 
the time interval, for the problem being considered, 
is divided. The .optimal ■ control satisfies the HJBE. It 
is assumed that the optimal control is unknown in the 
time interval being considered, but that a nominal 
control is available in the same time period. On 
application of this nominal control a nominal state 
trajectory, in the same time interval, is produced by 
equation (l). The nominal performance index value can 
be calciilated from equation (2). It is assumed that 
the optimal performance index expression is sufficiently 
well-behaved to allow a power series expansion in a 
small neighbourhood around the nominal state trajectory. 
In each iteration, the system equations are integrated 
in forward time, using the current nominal control. The 
accessory equations, which yield the coefficients of a 
linear or quadratic expansion of the performance index 
expression in the neighbourhood of the nominal state 
trajectory, are integrated in reverse time, thus 
yielding an improved control law because the performance 
index, which is a function of both state and control 
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varialsles, has been approximated within a small 
neighbourhood of the nominal state trajectory which 
was, in turn, produced by a nominal control. This 
improved control law is applied to the system equations 
producing a new and improved trajectory. By continued 
iterations, this procedure produces control functions 
that successively approximate the optimal control 
function. These successive approximation algorithms 
have been found to converge rapidly, for a large class 
of problems, to a locally minimal solution. Memory 
requirements of these algorithms are modest when 
compared with conventional dynamic programming. The 
sacrifice relative to conventional dynamic programming 
is, of course, global optimality. 

Another attempt to overcome the dimensionality 
barrier set up by conventional djniamic programming is 
the method of state increment dynamic programming. The 
essential difference between this technique and the 
conventional method lies in the choice of the time 
interval over which a given control is applied. In 
conventional dynamic programming, the sampling interval 
At is fixed, whereas in the state increment dynamic 
programming, the time increment 6t is the minimum time 
required for at least one of the state variables to 
change by one increment. Thus, the next state, for a 
given control, is known to lie within some small 
neighbourhood of the point at which the control is 
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applied. Phis procediore computes the optimal control in 
blocks -which may cover a long time interval but only a 
small distance along each state variable. 

The method of local variations is yet another 
example of a neighbourhood technique using a direct 
method of optimisation. Phis method has features of 
both the calculus of variations and djmamic programming. 
Chernous'Ko [56], and Erylov and . Chernous ' Ko [ 37 ] haye 
given a detailed description of the method and its 
applications. Prakasa Rao et al. [38] have applied this 
method to solve an optimal schediiling problem in hydro- 
thermal power systems. It has been found that this 
method may yield a local optimal solution and that 
its memory requirement is less than that required for 
incremental dynamic programming [38]. Banichuk et al.[ 39 ] 
have successfully applied the method of local variations 
to solve two point boundary value problems also. 

The salient features of this method are a 
discrete-time description of the system whose state 
equations are generally differential, an initial 
feasible trajectory as a starting point of the solution, 
and a systematic, iterative perturbation, in a specified 
small neighbo-urhood of the current feasible trajectory, 
such that the performance index value monotonically 
decreases from iteration to iteration till a satisfactory 
convergence to an optimal trajectory is obtained. At 
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every sta.ge of the algorithm, all the constraints of 
the problem are satisfied. In general, in each iteration, 
a sequence of mathematical programming problems heed to 
be solved. This aspect has been discussed in more detail 
in Chapter 2, Section 2.4- 

This method does not seem to have been applied 
widely for obtaining numerical solutions to optimal 
control iDroblems. There is, thus, a need to study this 
method in its application to optimal control problems 
with a view to gain computational experience and insight 
into the working of this method, as well as to see how 
it compares with the more widely used techniques of 
optimisation. 



CHAPTER 2 


THE METHOD OP LOCAL VARIATIOHS 

In tile previoiis chapter, the Method of Local Variations 
(MLV) has "been introduced as a neighbourhood technique using 
a direct method of optimisation. It has also been mentioned 
that this method exhibits features of both the Calculus 
of Variations and I) 3 mamic Programming, and m^y provide 
loca.lly optimal solutions to optimal control problems. In 
this chapter, a detailed description of the method will 
be given. This will be followed by its comparison with 
some Dynamic Programming techniques together with a list, 
with explanations, of the advantages and disadvantages of MLV, 

2.1 The Elementary Operation 

In order to facilitate understanding of the MLV 
algorithm, we will give a definition of an Elementary 
operation, ¥e will follow the approach of Krylov and 
Chernous' Ko [37] by formulating a typical optimal control 
problem to be solved by the MLV technique. 

Let us suppose that the process to be controlled 
is described by the 

State Equation; ^ f (x,u,t) (2,1-1) 

where x is the n-dimensional state vector, 

u is the m-dimensional control vector, 
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f is an n-dimensional vectop function 
wiiich is, in general, nonlinear and 
time-varying, with arguments as 
indicated, and 

t is the independent variable. 

The process is assumed to he described in the finite tine 
interval [o,T] which is considered to be fixed. 

It is required to minimise the performance index 


T 

J = / S (x,u,t) dt (2.1-2) 

0 

where g is a scalar fxmction of the arguments 

indicated in equation (2.1-2). 


subject to 

State Constraint ; G-(i) for o^ t^ T 

where &(t) is a variable closed region of the 

n-dimensional Cartesian space Oc 


and 


Contro3- Constraint; u(t)€ ir{t,x) for o^'t^T 


where 


U(t,x) is a variable closed region of the 


m-dimensional Cartesian space 


a 




under the 

Initial Condition; x(o) ^ G-(o) 
and 

Final Condition: x(T).€ G(T) 

let t^ and t 2 be two sufficiently close instants 
of time, and x^ and X 2 two sufficiently close state points 
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in the phase space, such that G-Ct^) and 

G-(t 2 ). Let us consider the problem of determining a 
control u(t) in the time interval [t^,t 23 , which translates 
the system given by equation (2,1-1) from co-ordinates 
(ti,Xi) to co-ordinates (t 2 ,X 2 ) and minimises the incremental 
performance index 


A J 





( 2 . 1 - 3 ) 


subject to the state and control constraints already 
stipulated being satisfied in the time interval [t^,t 2 ]. 
This problem will be referred to as an Elementary Problem , 
and the process of determining a solution to this problem 
as an Elementary Operation. 


Considering the pair of near-points (Ij.*— l^ 

^"^2^—2^* an elementary operation must provide: 

a) an answer as to whether there exists a control which 
translates the system given by equation (2.1-1) from 
co-ordinates (t^,x^) to co-ordinates (t 2 ,X 2 ) and satisfies 
the stipulated constraints ; 

b) if yes, an optimal control i^t) which solves the above 
variational problem, and the corresponding minimum incremental 
performance index value A jt Eor any pair - of sufficientl y 

Fo'f a-TvJ S>A.ffCo.e,in% 

^near-points .(t^f2i) elementary operation can 

be carried out approximately. Due to the proximity of (t^f2]_) 



and and tlie possi^bility of performing the elementary- 

operation approximately, the elementary problem becomes, 
in a large n-umber of cases, an algebraic problem only. In 
addition, the nearness of (t^,x^) and ('^ 2 *^ 2 ^ also be 
used for simplifying state equation (2. 1-1) by linearising 
it. The set G-(t) is usually a region of the n-dimensional 
state space changing continuously with time. As the ends 
of the trajectory, (t^,x^), (t 2 ,X 2 ), are sufficiently near 
and lie in this region, the state constraint set can be 
approximated by a linear inter polation between G'(t^) and 
G-Ct^) in [t^,t 2 ]. The set U(t,x) can be considered as 
approximately fixed, in carrying out the elementary operation at 

TJCt^x)^*^ [(t^+t2)/2, (xi+X2)/2] 

2.2 A Description of the ML7 Algorithm 

It is assumed that the elementary operation, 
satisfying the conditions form-ulated above, has been 
constructed for the problem under consideration. This 
elementary operation forms the basis of the Method of Local 
Variations algorithm; by the repeated, sequential application 
of this operation, it may be possible to obtain a locally 
optimal solution to the optimisation problem as described 
below. 

In order to apply the MLV algorithm, the state 
equation and the performance index expression, given by 
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eq.uations (2.1~l) and (2.1-2) respectively, are disoretised. 
This is done by dividing the time interval [o,T] into H 
equal subintervals to obtain a discrete-time representation 
of the system at time instants (k=so,l, . . . ,1T) , where 

i;=T/ 1I is sufficiently small. Next, an initial feasible 
tra.jectory x°(tj^); (k=o,l, , . . ,11) , is chosen such that the 
state and control constraints stipulated, including the 
boundary conditions on the state, are satisfied. The 
performance index value corresponding to this trajectory 
is calculated. 


This initial trajectory now becomes the current 
nominal trajectory at the zeroeth literation. Every state 
point on this trajectory is given a sequential perturbation 
of +h or -h according to the MIV procedure described below, 
giving rise to a sequence of elementary problems which may 
be solved by performing a sequence of elementary operations. 
After all the state points in the current nominal trajectory 
have been perturbed and the resulting elementary problems have 
been solved by elementary operations, one iteration is 
complete and a state trajectory which is better than the 
previous nominal trajectory, together vfith the corresponding 
control sequence and performance index value, become available 
for the next iteration. This newly generated state trajectory 
now becomes the current nominal trajectory and the MIV 
perturbation scheme is applied to it. This procedure is 
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repeated "until a specified stopping criterion is satisfied. 

Thus, at every iteration, a feasible trajectory, together 
with the corresponding control sequence and perfonnance index 
value, are available. 

The variation of the current nominal trajectory 
(k=o,l, . . . , 1 ?) , by the method of local variations 
consists of perturbing the state at each value of k by an 
n-dimensional pert"urbation vector h, and, if the perturbed 
state is admissible , solving the resulting elementary problem. 
This procedure is applied at time instants tj^; (k?=o,l, , . . ,N) , 
sequentially. The approach suggested in Krylov and Cherno"us ' 

Eo [37] for the perturbation scheme for elementary operations 
is the following; 

Each of the n components of x(t^)' ; (k=o,l, . . , ,K) , are 
successively perturbed in a particular order - say the first 
component, then the second and so on. Eor pertiirbation given 
to each component of the state at time instant t^^, the 
associated elementary problems are solved, before proceeding 
to perturb the next component of the state, at the same time 
instant. For each of the components x^; (j == l,...,n), its 
own perturbation step h.>o, is given. let the current 
nominal trajectory in the pth iteration, f (b:=o,l, , . . ,E) , 

be varied by the MLV technique, as described below, to obtain 
the nominal trajectory for the (p+l)th iteration, x^'^^(tj^) ; 
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(k=Oyl, . . . jjf) , Let it Le ass'uaed that the variation 
of the state values at the tine instant t£ is under 
consideration. Then variation of state values at tine instants 
t^; (^= 0 , 1 , . . , , 1 - 1 ) , have already heen performed . So the 
cui’rent nominal trajectory in the pth iteration at this 
sta e is; 

x^(tj^); (k=o,l, . . . , 1 - 1 ) , ^(tj^); .. ,N), where 

x^(tj^) ; (k=o,l, , , . ,1-1) , represents the new values 
assigned to the nominal state trajectory in the pth iteration 
hj the applico.tion of the MLV technique at tine instants 
t^; (k=o,l, . . . ,C-l) . Each coraponent of x^(tj^), starting 
from the first component and going upto the nth, is varied 
as follows, assuming, for simplicity, that all the components 
of the n-dinensional perturbation vector h are equal to h : 

A perturbation +h is given to the component under 
consideration. The following steps are then taken: 

a) A test is performed to see whether this perturbed value 
of the component under consideration of the state s^(\) is 
admissible? if admissible we go to step b) below, otherw3.se 
this step is repeated with a perturbation of -h to the same 
component? in the event of both perturbations leading to an 
inadmissible state value, we go to stop d) below? 

b) the optimisation problem of determining admissible 

control inputs in the time intervals [t 'V. ? 



20 


•wliich. take the system from the state ) at tine 

instant t^^ to the state x^(t|^j) at tine instant t^ via 
the perturbed state at time instant t^ , is then solved. 

If there does not exist a solution to this problem, the 
perturbation given to x^(t^ ) is considered a 'failure' and 
the above steps are repeated for pertiirbation in the sane 


component of x^(t^) of -h; if we meet with this 'failure ' 
for --h perturba,tion also, we go to step d) below; 


c) the increnenta.1 performance index values for the above 
state transfer in the time intervals [t^^ 

o.re calculated. If the sm of these two incremental perfor** 
nance index values is less than the corresponding sxm that 
exists before the above perturbation scheme for the state 
component under consideration, the perturbation is considered 
to be a 'success' and the perturbed state component value 
is assigned to the nominal trajectory at the time instant t^ 
and the whole perturbation scheme is then applied to the 
next component of the state at tine t^^ , 

d) If perturbations of both 4-h and -h meet with 'failures', 
the unperturbed value of the state component at time instant 

t|^ is retained and the above steps are repeated for perturbations 
of the next state component at the sane time instant t^ , 
till the state components have been subjected to the above 
operations . 
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After all the n conponents of the state ) 

have been varied by the above procedure, the algorithn 
moves on to the state point ) and varies the n 

coiaponents of that state point in the sane manner, In this 
vay, the local voxiations operation is performed successively 
for all n conponents of all the state points, upto time 
instant t^^. This completes the pth iteration, and the 
nominal state trajectory and control sequence, together with 
the modified perforjaance index value, are obtained for the 
(p+l)th iteration. The variation of the initial and final 
IJoints in the state trajectory differs in the facf that the 
elementary operations are carried out only " to the right” 
or '»to the left” respectively. Thus, each iteration 
generates a feasible state trajectory (alongwith the 
corresponding control , sequence and performance index value) 
which is better than (or at least as good as) the state 
trajectory in the previous iteration, in the sense of the 
performance index. These iterations, which successively 
reduce the perforaance index value, are performed till some 
convergence criterion is satisfied, thereby producing a 
locally optimal solution to the problem being considered. 

Pigure 2.2 explains the working of the local variations 
scheme for a scalar case. It shows a part of the nominal 
trajectory at any iteration and considers state values at 
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time instants (l-l), I and (l+l). A perturbation +h lias 
been given to the state value at the Ith time instant. 


XCHfiCK 



XCHSCK is the state value at the Ith time instant, perturbed 
by +h. 

11(1-1) is the nominal control input in the time interval [l-l,l]. 
U(l) is the nominal control input in the time interval [‘I,l4-l), 
DSLCST(l-l) is the incremental performance index value 

associated with the transfer of X(I-l) to X(l), 
DELCST(l) is the incremental performance index value associated 
with the transfer of X(I) to X(I+l). 

■ra'lIlJUS is the modified control input required to transfer 
X(l-l) to XOHECK. 

UPLUS is the modified control input required to transfer 
XCHECK to X(.I+1). 

DLCSTl is the modified incremental performance index value 
associated with the transfer of X(I-l) to XCHECK* 

I)1CST2 is the modified incremental performance index value 
associated with the transfer of XCHECK to X(I+l). 
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Tlie perturlDation given to X(I) poses an elementary 
problem. An elementary operation tries to solve this 
problem by finding out the following; 

a) whether the perturbed state XCHECK is admissible^ 

b) if yes, whether there exist control sequences UMIEUS 
and TJPLUS, capable of transferring X(I-l) to X(I+1) via 
XCHSCh; if so, whether IMINUS and UPIIJS are admissible j 
if yes, whether they are unique; if not, we are faced 
with a mathematical programming problem whose solution 
provides the optimal control inputs IMINUS and UPlUS; 

c) ass-'orning IJI‘aIiTUS and UPlUS exist, are admissible and 
unique, whether the snm of incremental performance index 
values (D1CST1+DICST2) associated with the transfer of 
X(l-l) to X(I+1) via XCHBCK, is less than the corresponding 
sum [DELCST(I-1)+DELGST(I)] associated with the transfer 

of X(I-l) to X(I+1) via X(I). 

If the answers to a), b) and c) above are all in the 
affirmative, the perturbation +h to X(l) is a 'success*. 
Otherwise, it is a 'failure* and a perturbation -h is now 
given to X(I). This situation can similarly be analysed as 
ab ove . 

2,3 Comparison of the Method of Local Variations with 
Dynamic Programming Techniques 

Bellman's Dynamic Program m ing [27,28,29] is a direct 
computational procedure for solving multi-stage decision 
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making problems by the repeated, sequential, stage-by-stage 
application of the principle of optimality. This technique 
provides a global optimal control law in a closed-loop form. 

It is applicable to both linear and nonlinear systems. As 
it is a direct search technique, any state or control 
constraints that may be stipulated in the problem, help to 
reduce the domain of search and ease the finding of a solution. 
The major limitation of dynamic programming is the storage 
requirement in digital computation. Bellman refers to this 
difficulty as the 'curse of dimensionality'. This feature 
is the penalty that must be paid for avoiding the TPBVP's 
that arise in the case of variational methods. Furthermore, 
as it is a systematic, exhaustive search procedure for 
multi-stage processes, the amount of computations may become 
prohibitively large for all but very low order systems. 

To avoid the 'curse of dimensionality', several 
modifications were done on the original dynamic programming 
algorithm. This resulted in the development of such techniques 
as Incremental Dynamic Programming [30 ], Differential Dynamic 
Programming [3l] and State Increment Dynamic Programming [3l], 
all of them being neighbourhood techniques using a direct 
method of optimisation. The Method of local Variations is 
yet another neighbourhood, technique utilising a direct method 
of optimisation and requiring significantly less storage 
capacity than the conventional dynamic programming technique. 



25 


Here also, any state or control constraints present in tlie 
problem help in' finding a sol'ution. The MI7 techniq.ue may 
■[irovide a local optimal control la'w and has been shown by 
Prakasa Rao et al.[38] to require a significantly less storage 
capacity than the incremental dynamic programming technique. 

Prakasa Rao et al.[38] have compared Incremental 
Dynamic Programming (IDP) with the Method of Local Variations 
in terms of computational requirements. Both IDP and HV 
are methods involving direct search in a neighbourhood of the 
nominal trajectory. In IDP, the principle of optimality 
is applied to obtain the best trajectory in- this nei^bourhood 
and the iterative procedure indicated in Section 1.2 enables 
a locally optimal trajectory to be obtained. On, the other 
hand, the MLV tries to obtain, in each iteration, only a 
better trajectory than the existing nominal trajectory, in 
the neighbourhood "under consideration. The incremental 
dynamic programming algorithm was suggested by Bernholtz 
and G-raham [30] in the context of optimal scheduling of 
hydro-thermal power generating systems. For the sake of 
comparison, Prakasa Rao et al.[38] considered a power system 
consisting of one hydro plant and one thermal plant. Let 
the interval of optimisation [o,T] be divided into H equal 
subintervals, each of duration At=l. 
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The hydro-system is described by the 


State Egtiation x(k+l) = x(k)+L(k)-u(k)~e(k) ; (k=l, , . , ,1T) 

(2,3-1) 

where x(k) is the storage capacity of the hydro 

reservoir at the beginning of the kth 
time period, 

L(k) is the inflow rate into the reservoir 
during the kth time period,. 
u(k) is the discharge rate from the reservoir 
d^iring the kth time period and 
e(k) is the evaporation rate from the reservoir 
during the kth time period. 

It is required to find the discharge rate u(k) such that the 
fuel cost in thermal power generation over the interval 
[o,T], given by 

J = S F (Frn(k)) (2.3-2), 

k=l 


is minimised subject to the 

Equality Constraint; Pj^(k)-Pjj(k)-Py(k)+Pj^Qgg(k)=o; (k=l, . . , ,K') 

(2.5-3) 


and the 

Inequality Constraints; 




,N) 
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inhere F(P^(k)) = a P^(k)+b Pj(k),* (k=l,...,N) (2,3-4), 

is the fuel cost of generating P^ thermal units during 
the kth time period, and ayh are constants, 

Pjj(k) is the load demand on the system in the kth 
time interval, 

P^(k) is the power generated in the hydro station in 
the kth time interval, 

P^(k) is the power generated at the thermal station 
during the kth time interval and 

PTr^c.„(k) is the power lost in transmission during the 
JuUoS 

kth time interval, 

The initial and final reservoir storages .x(l) and s^JJ+l) are 
specified . 

Prakasa Rao et al.[38] have considered a case in 
which the interval of optimisation for this one-hydro-one- 
thermal problem is divided into three equal subintervals 
as shown in Figure 2.3. 



x X 


Figure 2,3 
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The points a,lD,c,d represent the nominal state trajectory- 
in a particular iteration. The states e,f,g,h define the 
neighboiurhood considered in these methods. It has been 
sho-t/m that for IDP method, nine trajectories need to he 
compared to find the optimal trajectory in this neighbourhood, 
whereas for ML¥ technique, five trajectories are to be 
compared in the worst case, and three trajectories in the 
best case. The MIV technique not only requires a lesser 
number of trajectories for comparison than the IDP technique 
but also the nunber of main storage locations required by 
the MLV technique is significantly less than that required 
by the IDP technique. These results, obtained by Prakasa 
Rao et al.[38] for the above problem, have been tabulated in 
Tables 2.3(a) and (b); details can be found in Prakasa 
Rao et al. [38], 


Number of 
Intervals 

r 

IDP 


Q Trajectories to fi 
5 be compared I 

j_ 5 

Primary storage} No. of main 
requirement | calculations 

0 per iteration 

3 

9 

24 

9 

6 

243 

51 

24 

N 

3N-I 

3(3N-1) 

(5N-6) 


Table 2, 3 (a) 
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Number of 
Intervals 

T 

"r-“ 

MLV 



^Trajectories tS 

1 be compared | 

1 i 

Primary Storage > 
Requirement 

po. of main cal- 
jculations per 
iteration 

(jWorst ; 
lease I 

Best g 
case 0 

1 V/orst 
! case 

Best 

J case 

5 

5 

3 

15 

8 

•4 

6 

11 

6 

24 

20 

10 

N 

(2N-1) 

N 

3(N+2) 

4(N-1) 

2(N-1) 


Table 2.3(b) 

I’rom Tables 2.3(a) and (b) it can be concluded that 
for the problem under consideration, ML7 requires a significantly 
lesser number of storage locations, lesser number of trajectories 
for comparison, and fewer nuiaber of calculations per iteration 
than IDP. Thus, for the same nei^bourhood, ML7 seems to be 
preferable to IDP as a method of optimisation. When the 
initial trajectory is far away from the optimal one, MLV 
provides a faster transition from the current nominal trajectory 
to the next one, as compared to IDP, If the optimal trajectory 
lies within the neighbourhood being considered at any stage 
of the iterative procedures of IDP and ML7, IDP leads to the 
exact optimal trajectory, whereas MLV may provide only a 
better trajectory than the existing one. It is thus advisable 
to use MLV initially, till we are close to the optimal 
trajectory, and then switch over to IDP for convergence to 
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the optimal trajectory. Alternately, the disparity 
between the solutions given by IDP and MLV can be reduced 
by choosing progressively smaller perturbation steps in 
MIV as the optimum is approached, 

2 » 4 Some advantages and disadvantages of the 
M ethod of Local Variations 

listed below are some of the advantages and 

disadvantages of the MLV technique. 

Advantages ; 

1) The Method of local Variations has been found to be 
attractive due to its simplicity, efficiency and ease of 
implementation. ' 

2) As MIV is a neighbourhood technique using a direct method 
of optimisation, even before convergence to an optimal has 
been achieved, a feasible and practically satisfactory solution 
may be chosen. Also, any state or control constraints present 
in the problem under consideration help in narrowing down 

the domain of search and ease the finding of a solution to. 
that problem, 

3) Inequality constraints can be directly handled without 
using any 'artificial' techniques like the penalty function 
approach, 

4) This technique requires significantly less computer storage 
space than the conventional dynamic programming and IDP tech.-' . 
niques because only the solution for the current approximation 
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lias to be stored. As this approximation is improved, 
it can he ‘erased* hy replacing it by the new approximation. 

5) Chernous ' Ko [36] and Banichuk et al.[39] have applied 
the MLV algorithm to solve linear and nonlinear boundary 
value problems defined for simultaneous ordinary differential 
equations. These problems may arise in mathematical physics 
and optimal control theory. 

6 ) In the original algoritlm, the time interval of the problem 
being considered, is divided into N equal subintervals. It 
may be possible to gain computational advantage in the MIY 
technique by partitioning the tine interval into unequal 
subintervals. Some other modifications of the basic MLY 
algorithm to improve its efficiency have' been given in 
Chapters 3,4. 

Disadvantages 

1) The starting point in the solution to an optimal control 
problem by the MLY technique is an initial feasible trajectory. 
Generation of this trajectory may prove to be difficult, 
depending on the constraints on the problem being solved. 

Also, the solution to a problem may be sensitive to the 
initial feasible trajectory chosen as shown in Solution 1(b) 

to Problem 1 in Section 3.1. 

2) There is no guarantee that the MLY technique will provide 
an optimal solution to a problem, as has been shown in 
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Solution 1(a) to Problem 1 in Section 5.1. If at all it 
provides an optimal solution, the results would be generally 
valid only locally. 

5) In this technique, since the search is confined to a 
particular neighbourhood, we have to go from one neighbourhood 
to another in order to converge to a- locally optimal solution. 
If the step size chosen is very small at the initial stage 
of the solution, the number of iterations required for 
convergence to a local optimal will be very large, making 
the algorithm inefficient, 'This can be avoided by taking 
a crude step size in the initial stage of the solution and 
gradually refining it as we approach the optimum, ¥e can 
also, simultaneously, divide the problem time interval into 
gradually smaller subintervals and then apply the MIV . 
algorithm. 

4) In general, we may encounter a mathematical programming 
problem for each +h or -h perturbation. Referring to 
I'igure 2.2 we see that with a perturbation of +h to the 
state value at the Ith time instant, that is, Z(I), we are 
required to find control inputs 'UI4I1TUS and UPLTJS satisfying 
the control constraint, which will translate the state 
X(I-l) to the state X(I+l) via the perturbed state XCBECE 
which should satisfy the state constraint, such that the 
incremental performance index value for this path, given 
by (D1CST1+DLCST2) , is less than that for the original path 
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X(I-l) — ^X(I) — >X(l+l), given by [DELCST(I-1)+DBLCST(I) ] 
This involved sequence of nathematical programing problems 
in every iteration may be tine consuming in certain 
nr ob lens * 



CHAPTER .3 


C OlgUTATIOHAl RESULTS 


The algorithm of the Method *f Ii»ca.l Variations 
(MLV) has been applied to obtain nvimerical solutions 
to certain optimal control problems. As a result, some 
important featiires of the algorithm have emerged. 


3.1 A Linear. Constrained. Optimal Control Problem 
Problem 1 

State Equation ; x(k+l) = x(k) + u(k) j 

(k = 0,1, . . . ,H-1) (5.1-1) 

It is required to minimise the performance index 


J = x^d) 


+ 


N-1 



[u^(k)+x^(k)] 


(5.1-a) 


subject to the following constraints s 

State Constraint ; 0 ^ x(k) 1«5 r (k = 0,1, . , . ,11) 
Control Constraint ; -l4u(k) ^ 1 j (k = 0,1, ... ,1^-1) 
under the following boundary conditions; 

Initial Condition ; x(o) = 1.5 
Einal Condition ; 3:(N) is free. 


In this problem, ¥ = 100. 

Solution 1(a) 

Since the state equation and the performance 
index are given as discrete equations, they can be 
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directly utilised in the Method of local Variations. 

The first step to be taken is the generation of initial 
feasible state and control sequences. These are sequences 
which satisfy the boundary conditions as well as the 
stipiiLated state and control constraints. This problem 
of generating initial feasible state trajectories and 
the corresponding feasible control sequences may be, quite 
complicated, depending upon the boundary conditions and 
the state and control constraints. In fact, this is one 
of the disadvantages of direct methods of optimisation 
like the ELY which require a starting, feasible state 
trajectory and the corresponding control input. However, 
for the problem under consideration, an initial feasible 
state trajectory can be easily chosen. It was chosen 
to bei 

x(k) = 1.5 f (k = 0,1,...,H) , (3.1-3) 

The control sequence corresponding to this state trajectory 
is 

u(k) = 0 j (k = 0,1,...,N-1) (3.1-4) 

Prom equation (3.1-4) it is seen that the control 
sequence generated from the choice of the initial 
feasible state trajectory, is admissible. 

The initial state trajectory and control sequence 
were utilised by a computer program utilising the Method 
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of local Yariations to find the optimal solution. The 
main part of the program calculates the performance 
index value for the nominal state trajectory and control 
sequence and obtains an optimal control input and the 
corresponding state trajectory by repeatedly calling a 
subroutine B1E0PB.. This subroutine varies the existing 
feasible trajectory using the MLY technique and returns 
a feasible state trajectory and the corresponding 
control input which are better than the preceding ones 
in the sense of the perfoimance index. The state and 
control constraints are to be checked in the local 
variations procedure in the EIiEjZiPR subroutinei This 
is done by calling a subroutine CHKPS written for 
this purpose. The convergence criterion chosen for 
this problem was that the difference between the 
performance indices for two successive trajectories 
must be less than lO”^. The iterations were stopped 
when this criterion was satisfied. 

The results have been tabulated in Table 3*1 (a) 
and the state trajectory, the control sequence, and the 
performance index versus iteration number plots have 
been shown in Figs. 3*l(a), (b) and (c) respectively. 
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Ini tia,l value of performance index * ^=227.250 


SI. 

Ho. 

Perturbation 

step 

h 

Ho. of 
iterations 

K 

Min. value of 

performance 

index 

Execution 
time on 
IBM 7044 
(secs) 

1 

0.1 

17 

5.960 ■ 

5.217 

2 

0.05 

30 

5.262 

9.100 

3 

0.01 

150 

5.040 

25.70 


Table 3.1(a) 


The stR,te trajectory shown in Fig-. 3.1(a) seems to be 
of an anonclous character because, after dropping from 
its initial value of 1.5 to close to zero very quickly, 
it again shoots up to its starting value at the final 
time instant, though a better result could be obtained 
if the control, which, starting from -1.0 and quickly 
going to zero as shown in the results presented, had 
remained zero upto the final time instant. However, 
the algorithm is inca.pable, in this particular 
situation, of bringing about this solution. This 
behaviour brings out an important feat\ire of the ML¥ 
algoritlim. For the problem and the initial state 
trajectory under consideration, the convergence is in 
the sense of the algorithm only and a claim that the 
converged solution is locally optimal, cannot be made. 



RAJECTORY 
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This is due to the way in which the MLY technique tries 
to find a better feasible trajectory in a specified 
neighbourhood of the current nominal trajectory - the 
search for a better trajectory consists in perturbing the 
state value;, at a particular time instant leaving st8,te 
values at. all. other time instants xmperturb -^d, starting 
from the time instant corresponding to N * X ♦ 
sense of this perturbation scheme, the algorithm fails, 
in the present case, to obtain a better state trajectory 
than the one shown in Fig. 3*1 (a) in its neighbourhood 
def.ined by the perturbation step size h. This, of course, 
does not imply that there does not exist a better 
tra,jectory in this neighbourhood - there may exist a 
bettor neighbourhood trajectory,, but the algorithm may 
fa.il to find it. As a matter of fact, if the perturbation 
scheme of MLV is applied at N = l,^,*.*, 100 to the final 
state trajectory shown in Fig. 3*l(a)» all perturbations 
turn out to be failures and the trajectory of Fig.3‘l(a) 
is retained. This cannot, strictly, be called convergence 
to the optimal. This important feature of the algorithm 
should be borne in mind when the MIV technique is used 
for determining. optimal control. 

Solution 1(b) 

How, we attempt to obtain an optimal solution 
to Problem 1 starting with a different initial feasible 
state trajectory. The initial feasible state trajectory 
and the corresponding control input chosen, are shown 
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in Figures 3»l{d) and (e) respectively. The results 
obtained have been tabiaated in Table 3.1 (b) while 
Figures 3 •1(d), (e) and (f) show, respectively, the 
final state trajectory, control sequence, and the 
performance index plotted against iteration numbers. 


Initial value of 

perf ornance 

index = ^^ni 

= 57.420 

SI. 

I'lO. 

Perturbation 

step 

h 

No, of It~ 
erations 

K 

Min. value of 

performance- 

index 

Execution 
time on IBM 
7044 (secs) 

1 

0.1 

14 

4.076 

3.26T 

2 

0.05 

30 

3.770 

6.267 

3 

0.01 

141 

3.650 

27.40' 


Table 3.1 (Is) 

In this problem, the fina,l result is in accordance 
with our expectations and the solution obtained here 
is much better than Solution 1(a). The value of 
here is much less, than that in Solution 1(a). The 
values of here, for different values of h, are 

bettor than their counterparts in Solution 1(a). The 
number of iterations required here are either less 
(in two cases) or the same (in one case) as that in 
Solution 1(a). Finally, the execution times required 
here are less (in two cases) and more (in one case) 
than that in Solution 1(a). These results can be 



STATE VARIABLE VERSUS TIME PLOT 
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attributed to tlie fact that the choice of the' Initial 
feasible state trajectory and control sequence in 
Solution 1(b) were better than the earlier choice in 
Solution 1(a). Considering Solutions 1(a) andl(b), 
it can also be observed that the result of application 
of IUjY algorithm depends on the choice of the initial 
feasible state trajectory and control sequence. 

It is noticed that, for both Solutions 1(a) and 
1(b) to Problem 1, decrease in .step, size h results 
in rapid increase in the number of iterations required 
for convergence; consequently, the computation time also 
increas.es rapidly with decrease in h. 

3*2 A linear, Tx^o- State Variable, Unconstrained . 

O'otimal Control Problem 
Problem 2 

^1 ~ ^2 • 2“1 ) 

^2 ~ ”^2 ^ ( 3 • 2 - 2 ) 

It is required to minimise the performance index 
1 

J = / [x^ + + 0.005 u2] dt (3.2-3) 

subject to the 

Initial conditions ; Xj(o) = 0 


^2(0) = -1 
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•Tills prololeia was considered by Hsieh [ 41 ] and iias been 
solved by lasdon et al.[25] using the Conjugate Gradient 
Method. 

3 • 2.1 Solution with the original MLV Algorithm 

The first step, is to discretise the state 
equations given by ( 3 . 2-1) and (3.2-2) and the 
performance index given by (3,2-3). The time interval 
[ 0 , 1 ] is divided into 100 equal subintervals, each of 
duration At = 0 .01. Thus, the discrete tine approxi- 
mation to the problem is: 

St ate Equations 

X, (k+1) = 0.01 Xp(k) + X-, (k) ; (k=0,l., . . . , 99) 

(3.2-4) 

x^dc+l) = 0.99 X 2 (k) + 0.01 u(k) ; (3.2-5) 

(k = 0,1,. ..,99) 

It is required to minimise 

99 

J = 0.01 ^ [x,^(k) + x^(k)'+ 0.005 u^(k)] 

k=o 

( 3 . 2 - 6 ) 

subject to the 

Initial Conditions ; x^(o) = 0 

x^{o) = -1 

This discretised model was utilised and the 
MIV technique was applied to obtain a locally optimal- 
solution. The stopping criterion used was that the 
ratio of the difference between two consecutive 
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perfcrmance index values to the larger performance index 
value of the two, must he less than 10 . The results 

ohta.ined ha,ve been tabulated in Table 3»2(a) which also 
shows the results obtained by lasdon et al.[253‘. The 
initial feasible state traj ectory. and control sequence 
chosen, the final state trajectory and control sequence, 
and the performance index versus iteration number plots 
are given in Tigures 3 •2(a), (h) and (c) respectively.. 

Initial va.lue of performance index = 1*094 


SI. 

llo. 

Pertur- 

bation 

step 

h 

Fo. of 
Itera- 
tions 

K 

Min. value 
of perfor- 
mance index 

'^min 

Execution 
time on 

IBM 7044 
(secs) 

Optimum value 
of performance 
index- ilasdon 
et al.{25]) 

1 

. 0.01 

96 

0.1283 

12.00 


2 

0.005 

190 

0.1004 

19.00 

0.07 

3 

0.001 

lOOT 

0.0770 

84.00 



Table 3. 2 (a) 


Since lasdon et al.[25] have not given the computation 
time in their solution, we have not been able to make a 
comparison, in terms of computation time, of the MLV 
with their method for this problem. 

3*2.2 Solution with a modified version of the MLV 
Algorithm 

A modified version of the MLV algorithm was 
next used to solve the above problem. The motivation 
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for doing this was to see whether a faster convergence 
to a solution could he achieved. In the original MLY 
algorithm, if at any iteration, the kth point of 
current trajectory was perturbed by +h say, and this 
led to a ’success’, then the state value, control - 
sequence and the incremental performance index values 
corresponding to the perturbation, replaced the old 
state value, control sequence and incremental perfor- 
mance index values respectively. Then the algorithm 
went to the (k+l)th point and repeated the local 
variations procedure till the whole trajectory was 
swept. The case where a perturbation of -h was a 
'success' was similarly dealt with. In the modified 
version of the algorithm, if, at the kth point of the 
current trajectory a perturbation of +h is a 'success' 
then, instead of replacing the old state value, control 
sequence and incremental performance index values by the 
values corresponding to the perturbation and going to 
the (k+l)th point, a perturbation of +2h is now given 
to the kth point. If this is a 'success' also, a 
perturbation of +3h is given to the same point and so on 
till a maximum perturbation of +5h. Thus, within a 
limit of +5h', the state, control and incremental perfor- 
mance index values corresponding to the 'most successful' 
perturbation replace the old values. Then only does the 
algorithm go to the (k+l)th point. This procediire is 
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the modification. Another fact to be noted is that 
whereas h = 0.001 produced the smallest' value of 
and h = 0.01 the largest when the original algorithm 
was iised, vice versa is true when the modified 
algorithm was used, at least ijj the present case. The 
best result achieved by the modified algorithm 
using h = 0.01 is, however, not as good as that 
obtained by Lasdon et al.[25]. 

3*3 Linear. Two-State Variable. Constrained. Ontimal 
Oontrol Problems 

Problem 3 

The first problem of 'the above category 
chosen was; 

Sto.t^ Equations ; = x^ (3* 3-1 ) 

^2 ” “^2 ^ ( 3 « 3 - 2 ) ^ 

It is required to minimise the performance index 
1 

J = J* [x| + xj + 0.005 u^"] dt (3.3-3) 

subject to 

State Constraint ; ~ 0*5)^ + 0.5^ 0 

and the 

In itial Conditions; x^(o) = 0 


X2(o) = -1 
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This problem was originally considered by Jacobson 
and lele [ 42 ] using the G-eneralised Gradient Method 
wherein, the gradient of J is calculated with respect 
to a set of independent variables rather than with 
respect to control. In this method, the choice of 
independent variables is dictated by the constraints 
on the problem and could result in different combi- 
nations of state and control variables as independent 
variables along different parts of the trajectory. 

Here, the directions of search are determined using 
gradient projection and the conjugate gradient methods, 
Mehra a,nd Davis [43] have used some of the state 
variables, instead of the control variable, as the 
independent variables to solve the above problem. 

They have also indicated the possibility of choosing 
a mixed set of independent variables consisting of 
sta.te and control variables. The switchover from 
one set of independent variables to another occurs 
when the constraint boundary is hit. 

3 . 3 • 1 Solution with the original MIV Algorithm 

The constraint on X 2 (t) is a state-variable 
inequality constraint. In the application of MLV to 
this problem, the state variable ^2^^^ given 
perturbations. The perturbed values of x^Ct) and u(t) 
would solely depend on the perturbations given to 
x^Ct). An initial feasible X 2 (t) trajectory was chosen 
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as the starting point of the solution su^h that the 
corresponding z. (t) trajector.y and the control history 

, inIo,ll, ^ 

^{■b^/.are admissible. Discretisation was done as in 
Section 3.2.1, with N = 100 and At = 0.01* Then MIV 
was applied to obtain a locally optimal solution using 
the same stopping criterion as in Section 3«2.1, Due to 
the presence of a state constraint, the choice of an 
initial feasible x^Ct) trajectory becomes important. 

The initial feasible state trajectories of Section 3* 2,1 
happened to satisfy all the constraints and were, 
thus, chosen for this problem also. The results obtained 
by KI.V and Jtehra and Davis [43] have been tabulated in 
Table 3»3(s!') while the state trajectories, control 
sequences, and performance index versus iteration 
number plots have been shown in Figures 3*3(3'), (b) 
and ( c ) . 

Initial value of performance index = = 1. 093-68600 


SI. 

No. 

Perturba- 

tion 

step 

h 

No. of 
Itera- 
tions 

E 

Min. value 
of perfor- 
mance index 

^min 

Execution 
time on 
IBM 7044 
(secs) 

“Optimum value 
of performance 
index, (Mehra 
and Davis [ 43 ]) 

1 

0.01 

56 

0.30021131 

10.00 


2 

0.005 

161 

0.23267423 

19.00 

^::0.1T 

3 

0.001 

898 

0.21874910 

94.00 



Table 3.3(a) 
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Jacobson and lele [ 42 ] nsed a slack variable to trans- 
form the constrained problem into an unconstrained one. 
The MLV technique does not require such a transforma- 
tion and considers the inequality constraint as it is. 
This reduces the complexity of solution. 

Although Mehra and Davis [43] obtained a lower 
value of Jj^ ^-y, th.an the results achieved by MIY, their 
x^Ct) trajectory did not satisfy the initial condition. 
Using MIY, the best results were achieved with h = 0.001 
in 898 iterations, and all the constraints were satis- 
fied. In order to ease comparison with the results 
obtained by Mehra and Davis [43]» Figures 3 »3 (g)» (b.) 
and (i), giving the state trajectory, control sequence, 
and the performance index versus iteration number plots 
respectively, for h = 0.001, have been drawn.- 

3 . 3*2 Solution with a modification of the original MLY 
Algorithm 

Problem 3 has been solved with a slightly modified 
version of the original MDY technique. The modification 
consists of the following; 

a) An initial perturbation step h was chosen. 

b) This h was utilised to find the final state 
trajectory and control sequence, the performance 
index value and the number of iterations required, 
with the same initial feasible state trajectory and 
control sequence as in Section 3*3*1* 
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c) li was then halved and the final state trajectory 

and control sequence of h) above, became the initial 

feasible state trajectory and control sequence. The 
MI¥ algorithm was then used with this h and these 
in3.tial values to generate the corresponding final 
state trp.jectcry and control sequence. 

The procedure c) above was repeated till a stopping 
criterion, depending on the h- value used in Table 
3.3,(a) of Section 3«3*1 with which we are cempa: ring 
this modified method, wa,s satisfied. Then- the final 
state trajectories, control sequence, performance 
index value and the total number of iterations required 
were printed out. 

Double precision arithmetic was used in the 
computations and the results for the above modification 
to MIV and the original MIV, with h = 0.005 and h = 0.001, 
are presented below for comparison. 

Initial valvie of performance index = == 1.093687 

^double precision result) 

Original MLV Algorithm (with double precision arithmetic) 


Perturba- 

No, of 

Min. value of 

Execution time 

SI, tion 

It era- 

performance 

on 

No, step 

tions 

index 

IBM 7044 

h 

K 

'^min 

(secs) 

1 0.005 

161 

0.2321809 

42.00 

2 0.001 

898 

0.2187232 

216.0 


Table 3.3(13) 
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Comparing with. Table 3*3(^)f ior h = 0,0o5aii<i h*0»001, 
it is seen that although the ntimber of iteratiois required 
for convergence remains exactly the same and the value of 
'^min same, the execution times have increased 

considerably due to double precision arithmetic. A 
program listing, giving a solution to Problem 3 with 
doLible precision arithmetic and h. = 0.001, has been 
provided, 

MLV Algorithm modified as described above. 


sa. 

•Initial per- 

Pinal 

Total 

Min. value 

Execution 

No. 

turbation 

pertur- lo. of 

of per- 

time on 


step 

bation 

itera- 

formance 

IBM 7044 


^i 

step h-, tions 

index 

“^min 

(secs. ) 

1 " 

0.01 

0.005 

114 

0.2321922! 

33.00 

2 

0.016 

0.001 

213 

0.1817503 

65.00 


Table 3.3(c) 


Comparing with Tables 3. ,3 (a) and (b), for h *= 0.005 and 

h = 0.*001, it is observed that the modification has not 

only significantly reduced the number of iterations 

required for convergence as well as the execution time, 

but it has also reduced the value of = 0.001. 

min 

The value of J . for h = 0.005, however, remains almost 

mill 

the same. There could be some other way of choosing the 
step size from iteration to iteration which could further 
reduce J . as well as computation time. 
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A SOLUTION TO PRCBLZHB WITH DOUBLE PRECISION ARITHMETIC AND H=:J.GC 
THE INITIAL NOMINAL TRAJECTORIES AND CONTROL SEQUENCE ARE FED IN 

DOUBLE PRECISION XI ( lOA ) , X2 ( lOi ) i U ( 100 ) , DSLCSTI 100) 

DOUBLE PRECISION ?RIMDX,A,B 
XI (!)= ■. 

X2{1 j=-l. 

DO 1 1=1,50 
X2 ( I + 1 ) =X 2 { I ) 

xiii + j )=': . i*x2(i)+x: (I ) 

1 CONTINUE 

XI { 5Z ) =0. :.Vi *XC (51 ) + X:;. ( 51 ) 

X2(52)=-0.99 
DO 2 1=52,1- „ 

X2{I+:i )=K2( I )+C.Ci 

XI n +1 ) =0.0 i=<=x2 ( I )+ x:. c i ) 

V CONTINUE 
DO 6 1 = 1,1-' 

U{ I) = 100.>i'X2(I + C)-99.=f=X2{I) 

6 CONTINUE 

PRINT ir-D 

100 FORMAT! 15 X,=4=THe INITIAL NOMINAL TRAJECTORIES ARE*) 

PRINT 101, XI 
Idl FORMAT! 2K,1 .:dI3. 5) 

PRINT 1?^1,X2 
PRINT 200 

200 FORMAT ! 15 X, ’i'THi: INITIAL NOMINAL CONTROL SEQUENCE IS*) 

PRINT li-l,U 

CALCULATION OF PERFORMANCE INDEX 

PRINDX=': . "> 

DO 3 1=1, 

DELCST! I) =0.01*1X1! I )*Xi( I ) + X2! I )*X2! I ) ) + 0.00005*UI I )*UI 1) 
PRINDX=PRINDX+DELCST! I ) 

B CONTINUE 

PRINT l;.2,PRINDX 

102 FORMAT! iOX, *THi-; INITIAL VALUE OF P.l. =#,D15.7) 

K=1 

A=PRINDX 

10 CALL ELEOPR !X1,X2,U,DELCST) 

B=A 

A=O.L 

DO 4 1=1,100 
A=A+DELCST! I ) 

4 CONTINUE 

I F ! K-K / *5 1 ) 3 •: 2 , 3 0 1 , 30 2 

501 PRINT 102, <, A 





I;.3 F0RMAT(15X,*lTERATiaN NO. =*, 15, l4Xt*COST =4,015.7) 

302 B=(B-A)/B 

IF(B.LT.l.D-6)G0T05 
K= K+ i 
GOTOin 
C 

5 PRINT 104 

1U4 FORMAT (15X,4THF FINAL STATE TRAJECTORIES ARE GIVEN BY*) 
PRINT 101 ,X1 
PRINT I0l,X2 
PRINT 105 

l;5 F0RMAT(i5X,*THL FINAL CONTROL SECUENCE IS GIVEN BY*) 
PRINT ls,I,U 
PRINT 106, K, A 

106 F0RMAT(15X, *THE- FINAL NO. OF ITERATIONS =*,I5//14X, 

1 *THE FINAL VALUE OF P.I. =*,D15.7) 

STOP 

END 

SUBROUTINE ELEOPR ( Xi, X2 ,U ,OeLCST ) 

C THIS S/R PERFORMS THE ELEMENTARY OPERATION 

C 

DOUBLE PRECISION XI ( 101 ) , K2( 101 ) ,U ( XOO ) , OELCSTI 100 ) 

DOUBLE PRECISION H, SAVE ,XiCHK,X2 CHK,UMI NUS,UPLUS 

DOUBLE PRECISION DLCSTl ,DLCST2,A ,B,C,2 

00 I I=2,1C1 

H=0.001 

SAVE=X2(n 

J=l 

c 

2 X2CHK=SAVE+'H 
Z=I-1 

C=8.*(Z*S.^I-(:?.5)*{ Z*0.l>l-0.5)~».5 
IF( (X2CHK-C).GT.O.)GOT04 
UMINUS=100.*X2CHK-99.«X2( i-l ) 

DLCSTl=fl. 01*(Xl ( I-l )*Xl{ I-l)+X2( I-l )*X2( I-l) ) 

I *UNI NUS*UM I JUS 

IFd.EQ.lOl )GQT05 
C 

UPLUS=in3-*X2( I+l )-99.*X2CHK 
DLCST2=0-ei*<Xl(I)*X!.(I ) + X2CHK*X2CHK) 
i +0.00005*UPLUS*UPLUS 

c 

A=DLCSTl+DLCST2 
B=DELCST( I-l )+DELCST( I ) 

IFIA.GT.B )GQT04 

XI (I )=0.01*X2(I-1 )+Xl(I-l) 

X2(I )=X2CHK 
N=I-1 

UtN)=UMINUS 
U( I )=UPLUS 



DELCST(N)=DLCST1 
DELCSK I) =DLCST2 
GOTOl 

4 J=J+1 

IF( J.GT.2)G0T01 

H=-H 

G0T02 

5 A=DLCSTI. 

B=DELCST( I-J. ) 

IF (A.GT.B )G0T04 

XI {n=o.oi*x2(i-i j + 

X2(I )=X2CHK 

N=I~1 

UIM)=UMINUS 
DELCSK N) =A 
i CONTINUE 
RETURN 
END 
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3 * 3»-3 Solution with, the same modification as in 
Section 5.2.2 . 

Finally, Problem 5 was solved with the modification 
detailed in Section 3»3»2. Figures 3*3(d), (e) and (f) 
shovr tliG state trajectories, control sequences, and the 
ijerformancG index versus iteration number plots respectively. 
From the results tabulated in Table 3*3 (d), similar 
conclusions as in Section 3* 2.2 can be drawn. 


Initial value of performance index = ® 1.094 


SI. 

No. 

Perturbation 

step 

h 

No. of 
iterations 

K 

Min. value of 
performance 
index Jain ' 

Execution time 
on IBM 7044 
(secs) 

1 

0.01 

56 

C.5054 

14.00 

2 

0 . C;C5 

88 

0.3125 

22.00 

3 

C.COl 

102 

0.3470 

46.00 


Table 3»3(d) 


Problem 4 

The second problem in the category being considered, 
and the last problem chosen in this chapter was t 


# 

State Eauations: x^ = x^ 

(3.3-4) 

Xg = -Xg + u 

(3.3-5) 

It is required to minimise the performance 

. , 1 ■ -i" 

index 

J = J' [x| + Xp + 0,005 u*]dt 

(3-3^6) 
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subject to 

State Constraint ; “ ^*5)^ + ^*5 ^ 0 

and tlie 

Initial Conditions ; x^(o) = 0 

^^(o) = -1 

This problem too, was originally considered by Jacobson 
and lele [42] and has been solved by Mehra and Davis [43] 
using the Generalised Gradient Method,. 

3.3*4 Solution with the original MIV Algorithm 

The problem being considered is the same as 
Problera 3 in Section 3*3 with the exception that the 
state inequality constraint is now on x^(t) instead of 
Z2('fc). For application of the MLV algorithm the problem 
was put into the discrete form as earlier, and a 
particular control sequence u(k) ; (k = 0,1, . . . ,K-l) , 
was chosen such that the corresponding state trajectories 
were feasible. These intial feasible state trajectories 
were different from those in Problem 3» as shown in 
Fig. 3*3(j). 

Discretisation was done as in Section 3»2.1 with 
N = ICO and 4i t = 0. 01. The MLV algorithm was applied, 
with the same stopping criterion as in Section 3*2.1, to 
obtain a locally optimal solution. Double precision 
arithmetic was used for greater accuracy and the results 
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obtained have been tabulated in Table 3»3(e) together 
with the results obtained by Mehra and Davis [43]- 
Figures 3*3 (3)> (k) and (l) show the nature of the 
state trajectories, control sequences, and the perfor- 
mance index versus iteration number plots respectively. 

Initial value of performance index = *^ini ~ 454357 

Perturba,- No. of Min. value Execution Optimum 

SI. tion step Iterations of perfor- time on value of 
ITo. , „ mance index IBM 7044 performance 

^ ^ index [43] 


1 

■ C . Cl 

72 

0.5923687 

21.00 


2 

0.00 5 

189 

0.5146318 

48.00 

0.792942 

3 

C.CCl 

1C35 

0.5245792 

254.0- 



Table 3.3(e) 

Using MLV, the best results were obtained in 189 iterations 
with h = 0.005* Thus, in this application, the MIV has given 
a much better result than that obtained by Mehra and Davis 
[43]. Another aspect to be noted is that Figures 3.3(j) 
and (k) indicate that the final state trajectories and 
control sequences obtained with the help of MIY are quite 
different from those obtained by Mehra and Davis [43] » 
indicating that ML? has picked up a locally optimal solution 
different from that obtained by Mehra and Davis [43]* 
order to ease comparison with the results obtained by Mehra 
and Davis [43]? Figures 3.3 (p)f (q) and (r), showing the 
state trajectory, control sequence, and the performance index 
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versus iteration nunter plots respectively, for 
h = C.C05> have been drawn. 

3 • 3 • 5 Solution to Problem 4 with the MLY Algorithm 
modified as in Section 3.3.2 
Double precision arithmetic was used in the 
computations and the results for the above modification 
to MLY f with h = 0,005 aud h = C.OOl, are presented in 
lable 3»3(f) below for comparison with the results for 
the original MIV as shown in Table 3* 3(e). 


SI. 

TTo 

Initial 

Final 

Total No. 

Min. value 

Execution 

■ I'ertur- 

per tur- 

of itera~ 

of per- 

time on 


bation 

bation 

tions 

f ormanc e 

IBM 7044 


step 

step 

hf 

K 

ind ex 
“^min 

(secs) 

1 

C.Ol 

C.C05 

125 

0.5271251 

35.00 

2 

C.016 

0.001 

202 

0.5008761 

57.00 


Table 3.3(f) 


Comparing with Table 3*3(e)» for h = 0.C05 and h « O.OCl, 

it is seen that the modification has significantly reduced 

the number of iterations required for convergence and the 

execution times for both cases. The value of J . for 

min 

h = C.CCl has also been reduced while the value of 
for h = 0.C05 has increased slightly. 
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5 • 5 • 6 Solution with the ML7 Algorithm modified as in 
Section 5.2«2 

Finally , Problem 4 was solved with the above 
modification using double precision arithmetic. Figures 
3.3 (ei), (n) and (o) show the state trajectories, control 
sequences, and the performance index versus iteration 
number plots respectively, that were obtained. From the 
results tabulated in Table 3*3{g), ihe same conclusions 
as in Section 5*2.2 can be drawn. 

Initial value of performance index = ^=1*454357 


SI. 

ITo. 

Pertur- 

bation 

stop 

No. of 
Iterations 

K 

Min. value 
of perfor- 
mance index 

*^min 

Execution time 
on 

IBM 7044 
(secs) 

1 

C.Cl 

58 

0.6287411 

27.00 

22 

0.005 

74 

0.7870553 

37.00 

3 

0.001 

94 

0.8494794 

83.00 


Table 3.3(g) 

From the results obtained in Sections 3 . 2 . 2 , 3.3.3 and 
3 . 3 . 6 , the conclusion that can be drawn is that in 
terms of minimising the performance index, the ILV 
algorithm modified as in Section 3 * 2.2 does not give 
as good results as the original MLV algorithm. 
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In conclusion, it should he mentioned that In-this 
chapter, the MIV algorithm and some modifications on it, 
were applied to solve some optimal control problems. As 
a result, some computational experience regarding this 
technique, which was hitherto lacking, has been gained. • 
What is more, several important features of the MIV 
algorithm have emerged owing to its application to the 
problems solved. These traits of the algorithm can 
form a basis for its comparison with the more widely 
used techniques of optimisation. 



CHAPTER 4 


COECLUSIOHS AED SUGGEST lOHS POE EUTURE WORE 

In the previous chapters, the Method of Local 
Variations (MIV) has been studied and the original MIV 
algorithm, along with some modifications on it, have been 
applied to solve several optimal control problems. The 
computational experience gained in this study has helped 
us to recognise several of its important features which 
can provide a basis for its comparison with other better 
known techniques of optimisation. 

The Method of Local Variations has been described 
as a neighbourhood technique which utilizes a direct search 
procedure for optimisation. It is inherently simple and 
easily implementable , flexible and, thus, open to modifications, 
and efficient, in terms of memory requirement. Any state or 
control constraints stipulated in the problem ease the 
finding of a locally optimal solution,- These constraints 
can be directly handled by the MLV technique. However, as 
indicated in Solution 1(a) to Problem 1 in Section 3»1, 
this method does not always guarantee convergence to a local 
optimum. Another drawback is that the MLV technique requires 
an initial feasible trajectory to start with, whose generation 
may prove difficult depending on the constraints present in 
the problem. Indeed this problem is encountered in many 
direct methods of optimization. It has also been shown in 
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Solution 1(1d) to Problem 1 in Section 3fl that the result 
of application of this method depends on the initial 
feasible trajectory chosen. However, this is not a unique 
feature of this method; the solution obtained by most of 
the iterative optimization techniques depends dn the 
starting point for the iterations. 

Among the modifications attempted on the original 
MLT algorithm, the one reported in Section 3.3.2 gave 
promising results in terms of significant reduction in the 
number of iterations required for convergence and computation 
time. There may well be other modifications which may yield 
still better results than those obtained in Section 3.3.2. 

Por example, it would be interesting to see how the following 
MIV scheme works under the following conditions: 

a) The time interval of the problem under consideration 
is first divided into N equal subintervals and the MIV 
algorithm modified as in Section 3.3.2 applied to obtain 
a solution. 

b) With a starting point generated from the above solution 
(using an interpolation scheme), the same algorithm is 
applied to get another solution for a d is ere tiffed- time model 
of the continuous system defined over 2H equal sub intervals 
of the time interval of optimization. 

This procedure is repeated, dividing the time interval into 
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smaller and smaller subintervals till a satisfactory 
convergence to a locally optimal, solution is achieved-. 

We note in passing that there may be computational advantage 
in considering partitioning of the time interval of opti- 
mization into judiciously chosen unequal subintervals. 

In this thesis, only optimal control problems 
for linear systems have been solved by MLY, It is necessary 
to gain computational experience of M17 in its application 
to solve optimization problems defined for nonlinear systems* 

It would be worthwhile to investigate application of MIV 
to state equations linearized around the nominal trajectory 
and a simplified performance index defined only for the 
incremental state and control values'. Since MIV is essentia,lly 
a neighbourhood technique, we intuitively feel that such 
a scheme will work well, and there may be computational 
advantages to be gained. This linearizing idea can be 
found in differential dynamic programming technique [5l]» 

Prakasa Rao et al. [38] have suggested a way of 
getting over the difficulty which arises out of the uncertaintjr 
in the MLV algorithm to find a locally optimal soliition. 

When the initial nominal trajectory is far away from the 
optimal one, MIV can be applied to provide a faster transition 
from the current nominal trajectory to the next one. Once 
the current nominal trajectory enters a small neighbourhood 



60 


of the optinal one, we can switch over to the Incremental 
Dynamic Programming technique [30] which leads to the 
exact locally optinal trajectory. 

Finally, we suggest a modification regarding the 
apxolication of the MDV technique to solve multi-dimensional 
problems. The state trajectory, in this case, consists of 
a set of state-vector points at different instants tof time. 
According to the original algorithm, the first component 
of the state-vector, at a particular time instant is first 
perturbed by +h. If this is a ’success the corresponding 
substitutions in sta,te, control and incremental performance 
index values are made. If it is a ’failure that component 
is then perturbed by -h. If it is a 'success ’» the necessary 
siibsti tut ions are made, otherwise, the old state, control 
and incremental performance ihdex values are retained and 
the algorithm moves on to the second component of that state- 
vector, and so on, till all the components of that state- 
vector have been perturbed. Then the algorithm moves on 
to the state-vector at the next time instant and repeats 
the above process. This whole procedure is repeated till 
all the state-vector points, in the current nominal trajectory, 
have been subjected to the elementary operation. 

The suggested modification consists in ordering 
the above search. For this, the past history in terms of 
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• s’access-f allure ' patterns, of the system state has to 
he renemhered. If, say, the jth component of the state- 
vector at the ith tine instant and kth iteration had been 
successfully perturbed by -h, then it would be " logical” 
to assume that, as MIV is a neighbourhood technique, the 
chances that a perturbation of -h would be 'successful' 
for the same component of the state-vector at the sane 
tine instant, in the (k+l)th iteration, are much higher 
than those for a perturbation of +h to that component. This 
could also be true for the jth component of the state-vector 
at the (i+l)th tine instant and kth iteration. Thus, for 
this component of the state-vector at the (i+l)th time 
instant and kth iteration, and the ith time instant and 
(k+l)th iteration, we can directly give a perturbation of 
-h instead of first giving a perturbation of +h and then -h 
as in the original MIV algorithm. This should result in 
a significant reduction in execution time although this 
gain would be at the cost of an increase in memory require- 
ment needed to "remember" the past 'success-failure' 
patterns . 

From the above it is clear that many modifications 
of the basic MIV algorithm are possible, and there may be 
some which significantly reduce computational requirements. 
These are the possibilities that we see for extension of 
the work reported in this thesis. 



REFERENCES 


1. Atiiane , Status of Optimal Control Theory 

Applications for Eeterninistic Systems’, IEEE Trans., on 
Automatic Control, Vol, AC-11, No. 3, PP 580-596, July, 1966. 

2. Kirk, D.E., Optimal Control Theory 

An Introduction. Prentice Hall, Electrical 
Engineering series. 

3» Sage, A.P,, Optimum Systems Control . Prentice Hall, 

Electrical Engineering Series,,, 

4. Pierre, D.A., Optimisation Theory With Applications .- 

J ohn Wiley and Sons Inc . 

5. Bellman, R., and Kalaha, R., Quasilinearisation and Nonlinear 

Boundary Value Problems . Elsevier Press, 

New York, 1965. 

6., Bellman, R.E., Kagiwada, H.H., and Kalaha, R,B. , 'Quasilinear- 
isation, Boundary Value Problems and linear 
Programming', The Rand Corporation Memo RM- 
4284-PR, September, 1964. 

7. Bellman, R., Kalaba, R., and Sridhar, R. , 'Adaptive Control 

via Quasilinearisation and Differential 
Approximation* , The Rand Corporation, RM-3928-PR, 
November, 1963. 

8. Kalaba, R., 'On Nonlinear Biff er out: tal Equations, The 

Maximum Operation, and Monotone Convergence’, 
Journal of Mathematics and Mechanics, Vol,8, 

PP 519-574, 1959. 



9, McGill, R* , and Kenneth, P, , 'Solution of Variational 

Prohlens hy Means of a Generalised Mewton- 
Raphson Operator', AIAA Journal, Vol»2, 
pp 1761-1766, 1964. 

10* Kopp, R,E. , and McGill, R,, 'Several Trajectory Optimisation 

Techniques; Part I: Discussion*, A.V. Balakrishnan, 
and l.W. Feus tad t, eds. Computing Methods in 
Optimisation Prohleas , Academic Press, F.T., 
pp 65-89, 1964. 

11, Moyer, H.G. , and Pinlcham, G. , 'Several Trajectory Optimisation 

Techniques; Part II; ilpplication' , pp 91-105 
of reference [lO], 

12, Kumar, K.S.P,, and Sridhar, R.,'0n the Identification of 

Control Systems by the Quasilinearisation Method ' , 
IEEE Trans, on Automatic Control, Vol, AC-9, 

Fo.2, pp 151-154, April, 1964. 

13,. Detchmendy, D.M., and Sridhar, R. , 'On the Experimental 

Determination of the Dynamical Characteristics 
of Physical Systems,' Proc. of the Rational 
Electronics Conference, Vol..21, pp 575-580, 1965.' 

14# Sage, A,P. , and Eisenberg, B.R, , 'Experiments in Fonlinear 

and Fonstationary Systems Identification via 
Quasilinearisation and Differential Approximation' , 
Proc. of the Joint Automatic Control Conference, 
pp 522-530, 1965. 



15« Henrici, P., Discrete Variable Methods in Ordinary 


Differential Equations . John Wiley and Sons 
Inc., New York, 1962. 

16. Sylvester, R.J., and Meyer, P, , ’Two Point Botmdary Value 

Prohlens by Quasilineo,risation' , Journal 
Society Industrial and Applied Mathematics, 
Vol.13, pp 586-602, June, 1965. 

17. Sage, A.P, , and Burt, R.W, , 'Optimum Design and Error 

Analysis of Digital Integrators for Discrete 
System Simulation', American Pederation of 
Information Processing Societies, Proc» P.J.C.C 
V01.2T, Part I, pp 903-914, 1965. 

18. Sage, A.P., and ;3mith, S.L., 'Real-Time Digital Simulation 

for Systems Control', IEEE Proc,, Vol.54, I'To.l2 
pp 1802-1812, December, 1966. 

19. Bellman, R. , et al., 'Invariant Imbedding and Monline.o.r 

Piltering Theory', Rand Corp., EM-4374-PR, 
December, 1964. 

20. Detchmendy, D.M., and Sridhar, R., 'Seq.uential Estimation of 

State and Parameters in Noisy, Nonlinear 
Dynamical Systems', ASf!E Journal of Basic 
Engineering, Vol.88, Series D, No. 2, pp 362-366 
June, 1966, 

21. Sage, A.P., and Masters, G.¥. , 'On-line Estimation of States 

and Parameters for Discrete, Nonlinear Dsnaamic 
Systems', Proc. of the National Electronics 



Conference, Vol,22, pp 677’-682, 1966. 

22. Sage, A.P., and Masters, G.¥,,' Identification and 

Modeling of States and Parameters of Unclear Reactor 
Systems', IEEE Trans, on Unclear Science, pp 279-285, 
Eebrnary, 1967. 

23. Sage, A.P,, and Ellis, T.W., 'Seqnential Sntioptinal 

Adaptive Control of Uonlinear Systems -Proc, of the 
Ua,tional Electronics Conference, Vol.22, pp 692-697, 

1966. 

24. Ellis, T.¥. , and Sage, A.P,, 'Application of. a Method of on- 

line Comhined Estimation and Control,' Proc. Southwest 
IEEE Conference, .^ipril, 1968. 

25. lasdon, L.S., Mitter, S.K., and ¥aren, A.D., 'The Conjugate 

Gradient Method for Optimal Control Problems', IEEE 
Trans, on Automatic Control, Vol, AC-12, Uo.2, pp 132-138, 
April, 1967. 

26. Bryson, A.E., and Ho, Y.C, , Applied Optimal Control . ¥altham, 

Mass: Blaisdell, 1969. 

27. Bellman, R.E., and Dreyfus, S.E., Applied Dynamic Programming ., 

Princeton, U.J: Princeton University Press, 1962. 

28. Bellman, R.E., and Kalaba, R.E., Dynamic Programming and 

Modern Control Theory . Uew York: Academic Press, 1965. 

29. Bellman, R.E., Dynamic Programming. Princeton, U.J: 

Princeton University Press, 1957. 



50, Bernholtz, B., and (Jralian, I.J., 'Hydro-Themal Economic 

Scheduling ' , Part I-Solution by Incremental 
Dynamic Programming, Trans, A.I.E.E. ,• I960, 79, 
(III), pp 921-952. 

51. Jacobson, D.H., and Mayne, D,Q., Differential Dynamic 

Programming . American Elsevier Publishing Coy, 

Inc., Dew York, 1970. 

32. Larson, R.B., 'Dynamic Programming with Reduced Computational 

Requirements', IEEE Trans, on Automatic Control, 
Yol. AC-10, pp 155-145, 1965. 

35. Rosen, J.B., 'The Gradient Projection Method for Nonlinear 

Programming. Part I. Linear Constraints', J. SIAM 
(I960), pp 181-217. 

34. ' Rosen, J.B. , 'Iterative Solution of Nonlinear Optimal Control 

Problems', J. SIAM Control Series A (1966), 
pp 223-244. 

35. Rosen, J.B,, 'Optimal Control and Convex Programaing' , 

Nonlinear Pro graining. J. Abadie, ed. New York : 
John Wiley and Sons, Inc., 1967. 

36. Chernous'Ko, P.L., 'A Local Variation Method for Numerical 

Solution of Variational Problems', TJ.S.S.R. 
Computational Mathematics and Mathematical Physics, 
Vol.5, No. 4, PP 234-242, 1965. 



37. Krylov, l.A. , and Cliernous ’ Ko, F.I, , ‘Solution of Problens 

of Optimal Control by the Method of local Variations’, 
U.S.S.R* Computational Mathematics and Mathematical 
Physics, Vol.6, Uo,2, pp 12-31, 1966, 

38. Prcdcasa Rao, K.S,, Prabhu, S.S., a,nd Aggarwal, R.P.,., 

'Optimal Scheduling in Hydro-Thermal Power Systems 
by the Method of local Variations ' , International 
Conference on Systems and Control, Aug., 30-Sept .,1, 

1973. 

39. Banichuk, H.V., Petrov, V.M., and Chernous ' Ko, F.I. , 

'The Solution of Boundary Value Problems by the 
Method of Local Variations j F.S.S.R. Computational 
Mathematics and Mathematica.1 Physics, Vol.6, Bo. 6, 
pp 1-20, 1966. 

41., Leondes, C.T., Ed., Advances in Control Systems , Vol.2, 

1964-1966, New York; Academic. 

42., Jacobson, D.H., and Lele, M.M., 'A Transformation Techniq.ue 

for Solving Optimal Control Problems with a State 
Variable Inequality Constraint’, IEEE Trans, on 
Automatic Control, Vol.14, No, 5, PP 457-464, 

October, 1969. 

43. Mehra, R.K., and Davis, R.E., 'A Generalised Gradient Method 
for Optimal Control Problems with Inequality Cons- 
traints and Singular Arcs ' , IEEE Trans . on iiutomatic 
Control, V0I.17, No.l, pp 69-79, February, 1972. 



